Heat diffusion in a two-dimensional thermal fuse model 
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We present numerical studies of electrical breakdown in disordered materials using a two- 
dimensional thermal fuse model with heat diffusion. A conducting fuse is heated locally by a 
Joule heating term. Heat diffuses to neighbouring fuses by a diffusion term. When the temper- 
ature reaches a given threshold, the fuse breaks and turns into an insulator. The time dynamics 
is governed by the time scales related to the two terms, in the presence of quenched disorder in 
the conductances of the fuses. For the two limiting domains, when one time scale is much smaller 
than the other, we find that the global breakdown time t r follows t r ~ I 2 and t r ~ T 2 , where I 
is the applied current, and L is the system size. However, such power law does not apply in the 
intermediate domain where the competition between the two terms produces a subtle behaviour. 

PACS numbers: 62.20.M-, 46.50.+a, 07.05.Tp 



I. INTRODUCTION 

Motivated by important aspects such as failure predic- 
tion and material improvement, resistor network mod- 
els [H-Q have been intensively used for numerical stud- 
ies of mechanical and electrical breakdown phenomena 
in disordered media. The quasi-static fuse model [l|, al- 
though very simple, reproduces the basic damage regimes 
in breakdown phenomena [5|, [6[ . In the infinite disorder 
limit [7] the percolation regime is present. The fracture is 
totally disorder driven, and fuses burns out randomly un- 
til global breakdown is reached. If the disorder is small, 
there is little precursor damage until a single fracture is 
developing from one end of the system to the other. This 
is the nucleation regime. 

However, the quasi- static fuse model contains no real 
dynamics, and does not capture time dependency in cor- 
relations caused by the local currents. Dynamical effects 
has been included in the fuse model to study the elastic- 
ity problem [8[ . Sornette and Vanneste [9|, [l0| developed 
a model for electrical breakdown and plastic deforma- 
tion, which they referred to as the dynamic thermal fuse 
model. The temperature of a fuse was governed by a 
general Joule heating term, i b /g, and a heat loss term, 
— aT, which can be considered as a simplification of a full 
spatial diffusion description. This model has been exper- 
imentally realized by Lamaignere et al. [11], and later 
Mukherjee et al. [12|, [l3| by studying electrical break- 
down of carbon-polymer composites. This shows that 
the thermal fuse model is able to capture some of the 
phenomenology of breakdown in disordered media, like 
critical behaviour in the breakdown time. However, the 
model does not take into account the correlations due to 
heat diffusion between fuses. Thermal interaction with 
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neighbour fuses has previously been studied by Pennetta 
et al. 0, EH hi a biased percolation model. But this 
model implies instant thermalization of the fuses, and 
thereby neglects time dependent effects. 

With the thermal fuse model as a base we introduce 
spatial heat diffusion in the system. We study how the 
interplay of quenched disorder, current enhancement ef- 
fects and heat diffusion give rise to time dependent ef- 
fects which may seem counterintuitive with respect to 
the quasi-static fuse model and the biased percolation 
model. 



II. MODEL 

Our simulations are based on the thermal fuse model 
proposed by Sornette et al. 0, [To| . The model consists 
of a square lattice oriented at 45° with respect to the 
two boundaries opposite of each other, which act as bus- 
bars (see Fig. [1]). Each bond in the lattice is an electric 
fuse which behaves like an ohmic resistor when intact. A 
voltage is applied over the busbars which induces a total 
current / in the lattice. To each fuse j we assign a con- 
ductance gj from a power-law distribution p(g) oc g~ 1+l3 . 
The conductances are generated from gj = xf, where Xj 
is a uniformly distributed random number in the interval 
(0.5,1.5), and B = f3~ 1 . We call B the disorder param- 
eter. The values of the conductances are set at the be- 
ginning and never changed, corresponding to a quenched 
disorder. 

Our model differs from the thermal fuse model by Sor- 
nette et al. in the sense that the term —aT is replaced by 
a spatial diffusion term. We also use a fixed 6 = 2, which 
corresponds to the Joule heating effect. There is no loss 
of heat at the boundaries, so the system will always reach 
global breakdown for / > 0. We then arrive at the fol- 
lowing heat equation for our model in two dimensions in 
the continuum limit. 
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dT r 

-frfay't) = ^(x,y)+DV 2 T(x,y,t), (1) 

with boundary condition dT/dn = on the top and 
bottom boundary. Periodic boundary conditions are ap- 
plied in the horizontal direction. The heat capacity is 
C = 1 for every fuse, and D is the thermal diffusivity. 
The evolution in time is given by explicit Euler integra- 
tion, and adaptive timestep is used. When a fuse reaches 
its temperature threshold T r , chosen equal for all fuses, 
it breaks and irreversibly turns into an insulator. This is 
the way the temperature field reacts back on the current 
field. All the fuses start at equal temperature. When a 
fuse burns out, the current redistributes itself instanta- 
neously in the network. The network will then heat up 
until another fuse burns out. The total current is kept 
constant during the fracture process, with / = 1 for the 
results herein, if not specified otherwise. The fuses inter- 
act with the 6 nearest neighbours through heat diffusion 
(see Fig. [1]). Note that if more than one fuse reach the 




FIG. 1. A small network of L = 4. The black fuses are the 6 
heat-exchanging neighbours to fuse j. 

threshold T r within the same time step, those fuses are 
broken before the temperature field is updated. This en- 
sures that the currents are instantaneously redistributed. 
The current distribution is calculated by the conjugate 
gradient method. 

III. TWO COMPETING TIME SCALES 

A series of final fracture patterns for different values 
of D and B are shown in Fig. O The fracture which dis- 
connects the network in two pieces is outlined, and the 
temperature field is indicated by colors. B = 1 yields a 
uniform distribution, while B > 1 gives a broader dis- 
tribution with a tail towards large conductances. Note 
the lower cutoff at 0.5 B , and upper cutoff at 1.5 s in the 
distribution. 

The time dependence is governed by the two time 
scales To = gCT r /i 2 and t\ = £, 2 /D. The competition 
between these time scales generates different domains of 



behaviour. For To << t\ the domain is referred to as 
the large-current domain. If we consider D = 0, the 
large-current domain is present for all values of B. By 
comparing the fracture patterns for B = 1 and B = 5 
we see that for B = 5 a few large cracks appear, while 
for B = 1 the pattern is more diffuse, with single broken 
fuses. Due to the broader distribution for B = 5, the 
small conductances will heat up much faster than the 
large conductances, and the cracks will grow along path- 
ways of small conductances. This is contrary to what the 
quasi-static fuse model gives, where increased disorder re- 
sults in more single broken fuses, and a diffuse breaking 
pattern. This difference can be attributed to where the 
disorder is applied in the two models. In our model the 
disorder is in the conductances, while in the quasi-static 
fuse model the disorder is in the thresholds of the fuses. 

As D is increased, an intermediate domain appears. 
Since the diffusion smoothens temperature differences, 
the disorder in the temperature field decreases, and the 
effect of the redistribution of the local currents becomes 
more apparent. This is visible in form of a less diffuse 
damage pattern, and a more localized crack (see Fig. [2] 
for D = 0.1 and D = 0.5). 

Sufficiently large values of D yields To » ri, and the 
small-current domain appears. This is dominated by the 
redistribution of the local currents, and is characterized 
by a single crack disconnecting the system, with little 
precursor damage. Similar behaviour appears for low dis- 
order in the quasi-static fuse model. 

The three domains are visible when we plot the break- 
down time t r as a function of / on a log- log scale. This 
is shown in Fig. [3l For both the small- and large-current 
domain we find a good fit with the power law 

t r ~ r 2 . (2) 

To recover this, it is sufficient to show that the order 
of local breakdown events is independent of the value 
of /. For the large-current domain we can neglect the 
diffusion term in Eq. [TJ Then we realize that the first 
fuse to break will not change when / is changed, and that 
the local currents are proportional to the applied current, 
i oc /, between each breaking event. It follows that rest 
of the breaking sequence will follow in the same order 
when / is changed. A more complete derivation of Eq. [2] 
is given in Ref. [lOj. 

In the small-current domain the temperature differ- 
ences are small, and the value of the diffusion term in 
Eq. [1] will only be significant in a short period of time 
after each breaking event. Hence, it is reasonable to as- 
sume i oc I between each breaking event, and we obtain 
Eq. [2] for this domain as well. 

The presence of the three domains for different values 
of / is in agreement with the thermal fuse model by Sor- 
nette et al [i|,[lQ[, and experimental results. Lamaignere 
et al. [ll| fitted the time-to-failure data with Eq. [2] for 
large values of /, and t r ~ (I — / c ) -21 for small val- 
ues of /. Mukherjee et al. [13| found a good fit with 
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FIG. 2. (Color online) Temperature distribution with fracture patterns at the moment global failure is reached. The fracture 
which disconnects the network in two pieces is outlined. A system size of L — 100 is used, (a) B — 1 and (b) B — 5. A fixed 
seed is used for the random number generation so that the effect of the diffusion is clear. 



t r ~ (I 2 /Ic - I)" 1 - By letting I c — ^ 0, which is the case 
in our model due to no heat loss, these fitting laws are 
consistent with Eq. [2] In the intermediate domain there 
is a crossover from the small-current to the large-current 
domain. In this domain the order of local breakdown 
events depends strongly on the value of /, and the time 
dependence is more subtle. Fig. 2] shows the breakdown 
time versus D. We find that larger values of D result 
in longer £ r , since high temperatures are more effectively 
smoothened out. This result is qualitatively in agree- 
ment with the biased percolation model Q. However, 
a probabilistic approach to the disorder was used, while 
our model uses quenched disorder in the conductances of 
the fuses. This difference manifests itself in the impact 
D has on the percolation threshold p c . This is shown 



in the inset of Fig. |4j Diffusion causes less disorder in 
the temperature field, hence we find that p c decreases as 
D increases. This is contrary to what the biased perco- 
lation model gives, which approaches ordinary percola- 
tion (p c = 0.5) for vanishing disorder in the temperature 
field. We also find that the difference in t r between the 
two limiting domains for D = and D » in Fig. 2] is 
smaller for B — 1 than for B — 5. This difference reflects 
the time scale, n, at which temperature differences are 
smoothened out, and t\ decreases with decreasing disor- 
der. 
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FIG. 3. Log- log plot of the breakdown time t r as a function 
of /. The solid lines correspond to Eq. [21 The data points 
are averaged over 10 samples for L — 100. 
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FIG. 5. Log- log plot of the breakdown time t r as a function 
of L. The solid lines correspond to Eq. [3] with s = 2. 
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FIG. 4. Time to failure t r as a function of D with 1=1. 
Inset: percolation threshold p c as a function of D. The data 
points are averaged over 10 samples for L = 100. 

IV. SCALING ANALYSIS 

Based on analogy to percolation theory [l5| , we assume 
the following finite-size scaling relation 

t r ~L s . (3) 

Fig. [5] shows a log- log plot of t r as a function of L. 
The total electric resistance R is independent of L in a 
homogeneous two-dimensional system. We assume this 
for our model also, and we get i cx 1/L. Since i cx / in 
both the small- and large-current domains, it follows that 
t r ~ L 2 . Hence, for small values of L, the large-current 
domain appears (to << n). A crossover to the small- 
current domain (to >> n) is observed for sufficiently 
large values of L with D > 0. 

However, in the intermediate domain i cx I is not a 
valid assumption, and I can not be replaced by 1/L in 
the function for t r . 



V. ROUGHNESS OF THE FRACTURE 

We define the width of the fracture surface as the 
height-height fluctuations w = ((z 2 ) — (z) 2 ) 1 / 2 , where z 
is the height from the bottom of the network to the fuses 
that belong to the fracture surface, i.e the crack which 
disconnects the network. Studies of the quasi-static fuse 
model in two dimensions establish that the width scales 
as w ^ L c , with C 0.7 within 10% accuracy for different 
threshold distributions [l6|, [I?) . 

We generated between 1000 and 100 samples of sizes 
from L = 10 to L = 100, for various values of B and D. A 
log- log plot of w as a function of L is shown in Fig. [6] The 
slopes of the linear fits give the exponents ( listed in Tab. 
H The roughness exponent £ seems to be independent 
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FIG. 6. A log-log plot of w as a function of L. The solid lines 
are linear fits to the data, and gives the roughness exponents 
shown in Tab. |U The data for different parameter values are 
shifted vertically for clarity. 

of which time scale is dominating for B = 5, i.e large 
disorder. We obtain £ = 0.75 ± 0.1, which is in the range 
of the global roughness exponent observed in the quasi- 
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static fuse model. For B = 1 the roughness exponent 
is highly dependent on D, and is approaching unity for 
increasing D. It means that the fracture is guided by the 
structure of the network, and this is a trivial regime of 
behaviour. 



VI. DIVERGENCE OF THE RESISTANCE 

In the study by Lamaignere et al pj] they found that 
the total electrical resistance R follows the power law 



R~(t r -t)~ a , (4) 

in a critical region close to t r . A value of o^d ~ 0.65 
was obtained. A log- log plot of R versus (t r — t)/t r 
for our model is shown in Fig. [71 The obtained value 
is a = 0.28 ± 0.05. This is in agreement with the re- 
ported value for the thermal fuse model [lOj, and seems 
to be independent of the thermal diffusivity and disorder 
parameter. However, the critical region for power law 
scaling moves closer to t r as the small-current region is 
approached, and as the disorder in the conductances is 
decreased |18j . 



VII. CONCLUSION 

We have investigated a dynamic thermal fuse model 
with heat diffusion. The time dependence is governed 
by the two competing time scales To = gT r C/i 2 and 
T\ = £, 2 /D. The breakdown time follows t r ~ I~ 2 in 
both the small-current domain (to >> t\) and large- 
current domain (to << t\). This is in agreement with 
experiments on electrical breakdown of carbon-polymer 
composites pj], Ell- In the intermediate domain, com- 
petition between the two time scales produces a more 
complex behaviour. A characteristic feature of this do- 
main is that increasing the thermal diffusivity lengthens 
the lifetime t r of the system. 

Heat diffusion introduces new subtleties to the ther- 
mal fuse model, which still remain to be investigated. 



TABLE I. Values of £ for different thermal diffusivities and 
disorder parameters. 
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FIG. 7. Log- log plot of R as a function of (t r — t)/t r for 
different values of B and D. The data points are averaged over 
50 samples of L — 100. The solid lines have slopes between 
0.26 and 0.30, giving a « 0.28. For D = 10 we see that the 
critical region is very small, close to t r . 

However, the power law behaviour in the divergence of 
the resistance and in the roughness of the fracture has 
proven to be robust, and seems independent of the ther- 
mal diffusivity and the disorder in the conductances. 
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